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Systems of soft-core particles interacting via a two-scale potential are studied. The potential is 
responsible for peaks in the structure factor of the liquid state at two different but comparable 
length scales, and a similar bimodal structure is evident in the dispersion relation. Dynamical 
density functional theory in two dimensions is used to identify two novel states of this system, the 
crystal-liquid state, in which the majority of the particles are located on lattice sites but a minority 
remains free and so behaves like a liquid, and a 12-fold quasicrystalline state. Both are present even 
for deeply quenched liquids and are found in a regime in which the liquid is unstable with respect to 
modulations on the smaller scale only. As a result the system initially evolves towards a small scale 
crystal state; this state is not a minimum of the free energy, however, and so the system subsequently 
attempts to reorganize to generate the lower energy larger scale crystals. This dynamical process 
generates a disordered state with quasicrystalline domains, and takes place even when this large 
scale is linearly stable, i.e., it is a nonlinear process. With controlled initial conditions a perfect 
quasicrystal can form. The results are corroborated using Brownian dynamics simulations. 
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I. INTRODUCTION 

In hard condensed matter systems, the structure of the 
crystalline states that are formed is largely determined by 
the strength of the bonds between the atoms or molecules 
in the system, the dependence of the bonds on the ori¬ 
entation of the particles and the packing of the particles. 
In general, thermal fluctuations and entropy are less im¬ 
portant, unless one considers a system near the melting 
transition. In contrast, entropy and temperature can be 
all-important in determining the structure of soft matter 
systems. 

For polymers in solution, the interactions between 
pairs of chains depend on a delicate balance between 
energy and entropy [T]. When the solvent is good, the 
polymer chains form an open structure and interactions 
between pairs of polymers are largely repulsive and en- 
tropic in origin. On the other hand, when the solvent is 
less good, the polymer exhibits a tendency to collapse. 
In a good solvent, the strength of the repulsion depends 
on how branched the polymer is. As a result the form 
of the effective interaction between polymers can be tai¬ 
lored and controlled via the polymer architecture. In 
star-polymers, for example, the effective interaction po¬ 
tential is determined by the number of arms on each star 

m- 

The effective interaction between soft polymeric 
macromolecules is also soft. Since the centres of mass 
need not coincide with any particular monomer, the effec¬ 
tive interaction potential between the centres of mass can 
actually be finite for all values of the separation distance 
r between the centres. In this paper we discuss the struc¬ 
ture, phase behavior and dynamics of a two-dimensional 
(2D) model system of such soft core particles. 

The model that we study consists of soft particles 
which have a ‘core’ plus ‘corona’ (or shoulder) architec¬ 


ture. The particles interact via the following pair poten¬ 
tial: 

V{r) = (1) 

where R is the diameter of the cores of the particles, and 
Rs > R is the diameter of the corona (or shoulder) of the 
particles. In addition to the two length scales present 
in the potential, there are two energy scales: the energy 
penalty for a pair of particle cores to overlap is e(l + a) 
and the energy penalty for just the coronas to overlap is 
ea, where a is a dimensionless parameter that determines 
the shoulder repulsion strength. 

The particular form of the pair potential in Eq. 0 
arises from considering the effective interaction between 
the centres of mass of certain dendrimers or star- 
polymers. For dendrimers, this potential applies if the in¬ 
ner generations of monomers are of one kind (hydropho¬ 
bic, say), while the outer generations are of another kind 
(hydrophilic, say). Similarly, if a star polymer is made 
of diblock copolymers, then with a suitable choice of the 
block length ratio, the effective interaction potential be¬ 
tween the centres of mass is expected to be of the form 
in Eq. Q [HE]. 

In Fig. [l] we display the phase diagram for a system 
with temperature= landRs/R= 1.855. The fig¬ 
ure shows the (a,po) plane, where po = and (N) 

is the average number of particles in area . This phase 
diagram is determined using density functional theory 
(DFT), which is described below. Ref. [6] provides a 
brief account of some of the work elaborated here. At low 
densities po the particles form a liquid state. However, 
as the density increases, the particles overlap and then 
freeze to form one of two different crystalline states. The 
crystals are unusual: they are of the so-called ‘cluster- 
crystal’ variety [UlTHls], related to the fact that the par¬ 
ticles have a soft core. In the cluster-crystal multiple 
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FIG. 1: (Color online) Phase diagram in the (po,a) plane, 
where po is the average density. The system exhibits three 
phases: a liquid and two crystalline phases. The ‘crystal A’ 
is hexagonal, with a large lattice spacing, while ‘crystal B’, 
which is also hexagonal, has a smaller lattice spacing. The 
small-dashed blue line is the linear instability threshold (spin- 
odal) for the uniform liquid. The large-dashed green horizon¬ 
tal line is where the two principal peaks in the static structure 
factor S{k) have the same height, while the pink dotted line 
terminating in a circle is the locus where the two peaks in the 
dispersion relation uj{k) have equal height. The circle marks 
the point where the smaller k peak disappears. 


particles occupy each lattice site. When the parameter 
a = 0, the system reduces to particles interacting via a 
simple soft potential with one length scale and one energy 
scale. This is the generalized exponential model with ex¬ 
ponent n = 8, or GEM-8 model fluid [71414]. In 2D, at the 
temperatures relevant here, this model exhibits just one 
hexagonal crystal phase, with lattice spacing ^ R that 
is approximately constant with increasing density (note 
that at very low temperatures, a series of isostructural 
phase transitions is expected nanH]). Similarly, when 
a ^ 1 the contribution from the core of the potential 
becomes negligible and the fluid is again approximately 
a GEM-8 system, but now the particles have the larger 
diameter Rg and a stronger repulsion energy. Thus, for 
large a, the system forms a hexagonal crystal with lattice 
spacing ^ Rg. We henceforth refer to this larger lattice 
spacing crystal as the ‘crystal A’ phase, and the smaller 
lattice spacing crystal as the ‘crystal B’ phase. 

When the difference Rg — R is small, one can pass 
smoothly from one crystal phase to the other as a is 
varied (i.e., there is only one crystalline phase). How¬ 
ever, when the difference is larger, as in the case dis¬ 
played in Eig. crystal A and crystal B are two distinct 
phases separated by a phase transition at a ^ 0(1) IZO]. 
Indeed, many of the interesting novel properties of the 
present model described below all occur in the regime 
when a ^ 0(1), because they stem from the competition 
between the two different length scales, R and i?s, and 
the two different energy scales, (1 + a)e and ae. 

Two of the most striking properties of this system are: 


(i) When quenched to certain regions of the phase dia¬ 
gram, where a ^ 0(1), the system can sometimes freeze 
to form states with quasicrystalline order, (ii) On exam¬ 
ining in detail the crystal A phase, again when a ^ 0(1), 
we find that there is a high proportion of mobile parti¬ 
cles in the system, which is why we refer to this phase as 
a ‘crystal-liquid’. It is crystalline, because the majority 
of the particles in the system are frozen onto a regu¬ 
lar hexagonal lattice. However, a minority are ‘liquid’, 
in the sense that they move throughout the system. We 
find that the proportion of mobile particles in the system 
can be as high as 7%. 

The quasicrystals (QCs) formed by the present sys¬ 
tem are a local equilibrium state of the system, i.e., they 
are not the global minimum free energy state. They are 
found at state points in the portion of the phase dia¬ 
gram where the thermodynamic equilibrium phase is the 
crystal-liquid A state. The QCs are formed by a particu¬ 
lar dynamic mechanism when the system is quenched to 
certain regions of the phase diagram. 

In order to find parameter values at which QCs might 
be favored, we have invoked understanding developed 
from analyzing mode interactions in the Earaday wave 
experiment, in which a tray of liquid is subjected to 
vertical vibrations of sufficient amplitude that standing 
waves form on the surface. This system exhibits quasi¬ 
patterns (the fluid dynamical analogue of quasicrystals), 
as discovered in the early 1990 ’s umin], and two dif¬ 
ferent mechanisms for their formation have been iden¬ 
tified (see HZ] for a more detailed discussion). Briefly, 
patterns with Q-fold symmetry are expressed as sums of 
modes with Q wavevectors spaced at equal angles, and 
weakly nonlinear theory is used to compute how waves 
with different orientations affect each other. One mecha¬ 
nism relies on strong self-coupling to downplay the effect 
of waves with different orientations [T8142Q] , so permitting 
8, 10, 12, 14, 16, 18, 20-fold or higher quasipatterns [T7| . 
The second mechanism invokes nonlinear coupling be¬ 
tween the primary waves with secondary weakly damped 
(or weakly excited) waves, such that primary waves with 
wavevectors separated by a certain angle determined by 
the ratio of the primary to secondary wavenumber, are fa¬ 
vored [laiiiHizi- We invoke here this second mechanism, 
as done in [28l |29] , and select the length scale ratio for 
our investigation to be Rg/R = 1.855 (Sec. HI) in order 
that the ratio of the primary to secondary wavenumbers 
is 2cos(15°) = 1.932, so favoring dodecagonal quasicrys¬ 
tals. 


In fact the mechanism for QC formation that we ac¬ 
tually observe differs from either of the two mechanisms 
described above. The QCs form when the uniform liq¬ 
uid is linearly unstable against density fluctuations with 
a small wavelength that is close to that of the lattice 
spacing of the crystal B phase but stable with respect 
to wavelengths comparable to that of crystal A. Thus, 
in the initial stages after a quench the system appears 
to be forming the crystal B phase. However, the mini¬ 
mum free energy structure is actually the larger lattice 
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spacing crystal A phase. In the subsequent nonlinear 
evolution, the system seeks to form this larger lattice¬ 
spacing phase. However, being already patterned with 
the shorter length scale from the early stage linear dy¬ 
namics, the system cannot always form a perfect crystal 
A and often forms a state with a mixture of both the 
short and long length scales that sometimes turns out to 
have quasicrystalline ordering, i.e., the Fourier transform 
of the density distribution reveals the presence of 12-fold 
ordering. As one might expect from such a mechanism, 
the structure that is formed contains defects. However, 
by carefully controlling the wave numbers of the density 
modulations prior to the quench, the system can be in¬ 
duced to form a ‘perfect’ quasicrystal. 

Understanding the mechanisms by which soft matter 
QCs can form is becoming increasingly important. The 
possibility of designing soft-matter quasicrystals that 
self-assemble has generated considerable interest at a fun¬ 
damental level, leading to a burst of experimental and 
theoretical activity [30]-[34]. Self-assembled soft-matter 
quasicrystals are of interest for a number of reasons, not 
least because they promise to provide a route to manufac¬ 
turing materials and coatings with novel optical or elec¬ 
tronic properties arising as a consequence of their high 
degree of rotation symmetry |35ti37] . 

Although our focus is on polymeric soft matter QCs, 
we should also mention that there are other (colloidal) 
soft matter systems that form QCs [38ti43] . These also 
have pair potentials involving more than one length scale, 
but owing to the particles having a hard core, the local 
structure of these materials differs from that described 
below, as does the resulting phase behavior. 

This paper is structured as follows: In Sec. |H| we de¬ 
scribe the DFT and dynamical DFT that we use to deter¬ 
mine the structures formed by the system. In Sec. |HI| we 
discuss the properties of the uniform liquid state, present¬ 
ing results for the radial distribution function g{r), the 
static structure factor S{k) and the dispersion relation 
uj{k). We then present results relating to the solid states 
that are formed, focusing on the crystal-liquid state in 
Sec. rV and on QC formation in Sec. |V| This section in¬ 
cludes a discussion of our numerical results and their rela¬ 
tion to other mechanisms of QC formation from the liter¬ 
ature that are relevant to soft-matter systems. The paper 
concludes in Sec. VI with a few concluding remarks. 


II. THEORY FOR THE SYSTEM 

We use DFT [44ll47] to determine the structure, ther¬ 
modynamics and phase behavior of the system. To de¬ 
scribe the dynamics of the system when it is out of 
equilibrium, we use dynamical density functional theory 
(DDFT) j48ll5T] . The thermodynamic grand potential of 
the system is a functional of the one-body density distri¬ 
bution p(r) of the particles: 

^[P(r)] = F{pi^)] + j dr/>(r)($(r) - /u), (2) 


where /i is the chemical potential, $(r) is the external 
potential and F[p\ is the intrinsic Helmholtz free energy 
of the system, which is composed of two contributions: 

F[p{r)] = ksT J dr/9(r) [ln(/9(r)A2) - l] + [/>(r)]. 

(3) 

The first term is the ideal-gas contribution, with ks the 
Boltzmann constant, T the temperature and A the ther¬ 
mal de Broglie wavelength. The second term, is the 
excess (beyond ideal gas) portion describing the contri¬ 
bution to the free energy stemming from the interactions 
among the particles. The equilibrium density profile of 
the system at a given state point (/i, T) is that which 
minimizes U[p], i.e., which satisfies the equation 

,4) 

5p(y) 

For systems of soft-core particles such as those we con¬ 
sider here, the following rather simple mean-field approx¬ 
imation is remarkably accurate [Ill52l|53]: 


FexW)] = ^ J dr J drV(r)y(|r-r'|)/>(r'). (5) 

This functional generates the following simple random- 
phase approximation (RPA) for the pair direct correla¬ 
tion function: 


c(2)(|r _ F\) = —/3 


S^Fexlpjr)] 

5p{r)Sp{r') 


-PVi\r-r'\), (6) 


where (3 = {kBT)~^. 

To calculate the density profile of the system at a given 
state point we discretize the density profile on 

a square Cartesian grid and use fast Fourier transforms 
to evaluate the convolution integrals in Fex[p]- We em¬ 
ploy standard Picard iteration [54] to solve the Euler- 
Lagrange equation obtained from Eqs. (§ @, 

ln[p(r)A2] - (r) + /?$(r) -pp = 0, (7) 


where 




SFex[p{r)] 

Sp{r) 


(8) 


is the one-body direct correlation function. Eor the RPA 
functional in Eq. this gives 

= — J dr';dU(|r — r'l)p(r'). (9) 

Equation Q can be rearranged to obtain the following 
expression for the density profile, 

p{r) = po exp (-^$(r) + (r) - [po]) , (10) 

where [po] denotes the value of when Eq. is 
evaluated for the uniform density profile p(r) = po. We 
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note the result po = A“^exp(—;5/i + c*^^^[po]) showing 
that the average density in the system po is determined 
by the chemical potential p or vice versa. 

Picard iteration of Eq. ( p!Q| ) corresponds to substituting 
the density profile at step j, p^-^^(r), into the right side 
of Eq. (10) to obtain To stabilize the iteration 


process these two density profiles at step j are mixed, 

p(j+i)(r) = a/)^£(r) + (1 - a)p^^'>{r), (11) 

to obtain a new approximation for the density at step 
j + 1. This equation is then iterated until convergence 
is achieved. The value of the mixing parameter a varies, 
depending on the state point and the type of density 
profile to be calculated, but typically is in the range 
0.001 <a<0.1. 

In the present study, we generate the initial guess for 
the density profile in several different ways. One choice 
is to use the density profile obtained from solving at a 
different state point. Another way is to start with the 
density profile p(r) = po + ^(r), where ^(r) is a small 
amplitude randomly fluctuating field. When the uniform 
fluid is linearly unstable (see below in Sec. IIC), this 
initial guess may converge to the density profile of the 
crystal. However, this method often results in density 
profiles containing defects. The chance of these forming is 
much less in smaller systems and so the density profile for 
a larger portion of a perfect crystal needs to be built up 
from the density profile obtained from a smaller system. 

With this procedure the average density in the system 


P = 


^ J drp{r) 


( 12 ) 


equals po only when the system is in the uniform liq¬ 
uid state. Eor the crystal, p 7 ^ po- This is because 
by iterating (10), we actually select the value of the 
chemical potential p. The density po is the density of 
the uniform liquid for this value of p, and p is there¬ 
fore the density of the crystal corresponding to this p 
value. In calculations where we wish to specify the av¬ 
erage density to be po, we add an additional step to the 
Picard iteration, where at each step j, after the mixing 
step given by Eq. we renormalize the density pro¬ 
file, whereby we replace with /p*^-^+^^(r), where 

f = Po/[^ J drp^-^+^^(r)], cf. Eq. ( p^ . In all our dis¬ 
cussions below, we do not distinguish between p and po- 
We use Po to denote the average density in all phases, 
but it should be borne in mind that when this refers to 
a crystal phase, we mean the average density as defined 
in Eq. (12). 


As presented, Picard iteration is simply a numerical 
algorithm for solving the Euler-Lagrange equation and 
therefore for finding density profiles which minimize the 
free energy. However, as shown in Sec. |V| Picard itera¬ 
tion generates a series of density profiles that are often 
a fairly good approximation t o the real dynamics as de¬ 
termined by DDET (see Sec. HB), i.e., the index j can 
be thought of as if it were proportional to the time t. 



FIG. 2 : (Color online) Top: density profile in the (x/R,y/R) 
plane near an interface between coexisting crystal A and crys¬ 
tal B phases, for a = 0.75. Because the lattice spacings of the 
two crystal structures are not commensurate defects along the 
interface are necessarily present. Bottom: density profile be¬ 
tween the uniform liquid and the [ 1 , 1 ] interface of the crystal 
B phase, for a = 0.5. 


In these cases we have used the fictitious dynamics gen¬ 
erated by Picard iteration in place of the slower DDET 
to survey the behavior at different points in the phase 
diagram. 

In Eig. we display typical density profiles obtained 
from DET when the temperature ksTje = 1 and Rg/R = 
1.855. The phase diagram for this system is displayed in 
Fig. [3 In the left panel of Eig. we display the density 
profile at an interface between coexisting crystal A and 
crystal B phases, for a = 0.75. Since the lattice spacings 
of the two crystal structures are incommensurate, defects 
are necessarily present along the interface. In the right 
panel of Eig. we display the density profile across an 
interface between the uniform density liquid on the right 
and the small lattice spacing crystal B phase on the left, 
for a = 0.5. 

In the phase diagram displayed in Eig. the shaded 
(red online) regions indicate coexistence between two dif¬ 
ferent phases. The boundaries of these regions corre¬ 
spond to the densities of the two phases at coexistence; 
recall that for two phases to coexist the chemical poten¬ 
tial /i, the pressure p = and the temperature T 

must be equal in the two phases. There is also a triple 
point, where all three phases coexist. Note that in order 
to determine the minimum grand potential Cl one must 
also minimize with respect to the computational domain 
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size L. Actually, for hexagonal crystals one should calcu¬ 
late the density on a domain of size x L. However, on 
comparing results on such domains with those obtained 
on a square domain, we have confirmed that as long as 
the (square) domain is sufficiently large that it contains 
many unit cells, the slight strain energy contribution to 
the free energy is negligible. Most of the calculations 
presented here are for a system of size L = 25.6i? with 
periodic boundary conditions, where the finite size effects 
for the regular crystal structures are negligible. For QC 
structures, there are particular domain sizes (e.g., 8, 30 
and 112 times the smaller scale) that allow for accurate 
approximation to 12-fold QC structure El- Our results 
are for domains of size 30. We remark further on this 
point in Sec. [V| 


A. Liquid structure factor 


To characterize the structure of the liquid state, two 
quantities, the real space radial distribution function g{r) 
and the reciprocal space static structure factor S(k)^ are 
very useful m- The static structure factor in the liq¬ 
uid phase is given by the relation S{k) = [1 — pQc{k)]~^^ 
where c{k) is the Fourier transform of Thus the 

RPA approximation for the structure factor in the uni¬ 
form liquid phase is 


S{k) =-.-. 

l + p^pV{k) 


(13) 


To determine the radial distribution function g{r) one 
can insert the simple RPA approximation ^ into the 
Ornstein-Zernike equation [47]. However, we choose in¬ 
stead to calculate g{r) using the Percus test particle 
method gZlES]. This gives a more accurate approxima¬ 
tion for g{r) and also illustrates better the true accuracy 
of the DFT that we use. The test particle method corre¬ 
sponds to fixing one of the particles at the origin, so that 
4>(r) = V{r) in Eq. ([^, and then calculating the density 
distribution of the remaining particles in the presence of 
this fixed particle. The radial distribution function is 
obtained from the resulting density profile via Eq. 
g{r) = p{r)/po. In Ref. [55| results from this RPA-test- 
particle theory were compared with the more sophisti¬ 
cated hyper-netted-chain (HNC) theory for a very simi¬ 
lar 2D soft-core system. The agreement between the two 
is rather good, which gives us confidence that the simple 
RPA DET is acc urat e. We present typical results for g{r) 
and S{k) in Sec. HI below. 


B. Dynamics: time evolution of the density 

In addition to the equilibrium fluid structure, we also 
determine the non-equilibrium fluid dynamics. Since we 
consider soft polymeric ‘blobs’ in solution, an appropri¬ 
ate approximation is to assume that the centres of mass 


of the particles move via Brownian motion, i.e., via over¬ 
damped stochastic equations of motion: 

ri = -rvc/(r^,t) + rx(t), (14) 

where i = 1,.., A" is an index that labels all the different 
particles in the system, whose set of position coordinates 
we denote by = {ri, r 2 , • • • , rjv}- The mobility coeffi¬ 
cient F = /3D, where D is the diffusion coefficient, while 
X(t) denotes the random force on the particles due to the 
solvent thermal motion. We assume in the standard way 
that X(t) is a Gaussian random variable j48ll5T] . The 
potential energy of the system is 

N N 

= y]$(ri) + y]y]l/(|ri -rjl). (15) 

i=l j>i i=l 


For a system of interacting particles with equations of 
motion given by ( [l4| ), we can use DDFT |48]-[5T] to de¬ 
termine the time evolution of the fluid non-equilibrium 
density distribution p(r, t). In DDFT the dynamics is 
governed by 


dp{r, t) 
dt 


= rv- 


p{r,t)V 


6Q.[p{r,t)Y 
5p{r,t) J ’ 


(16) 


a result that follows on making the approximation that 
the non-equilibrium fluid two-point density correlation 
function is the same as that in the equilibrium fluid with 
the same one-body density distribution. Equation (16) 
is thus an approximation 


but for soft-core fluids, 
previous good agreement with the results from Brownian 
dynamics (BD) computer simulations (i.e., from solving 
repeatedly Eqs. and then averaging over the differ¬ 
ent realizations of the noise), gives us confidence that 
Eq. (16) provides a good approximation to the exact dy- 


C. Dispersion relation 


An important quantity for understanding the behav¬ 
ior of the system is the dispersion relation, uj{k). This 
relation determines the rate at which density fluctu¬ 
ations in the uniform liquid grow {uj > 0) or decay 
{uj < 0) over time. Consider a uniform liquid with den¬ 
sity po: with a superposed small amplitude perturbation 
p(r,t) = p(r, t) — pQ. Equation (16) shows that the per¬ 
turbation evolves according to 


^^-Cp + Oip), 


(17) 


where C = is a linear operator and G 

denotes a convolution, i.e., c^^^(S)p = J (ir'c^^^(r—r')p(r'). 
To obtain this result one must make a functional Taylor 
expansion of Fex [p] uni E3 EE] • Linearising Eq. ( [T7| ) and 
decomposing p into a sum of different Fourier modes, 


p(r,f)=^Pke*'^-+“«*, fc=|k|, (18) 

k 











leads to the dispersion relation ISniESlESl 
w(fc) = -Dk^[l - poc{k)]. 


(19) 
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On combining this result with the RPA approximation 
we obtain = —Dk‘^[l-\-poPV{k)]^ a result closely 
connected to the structure factor S{k) defined in Eq. (13). 
This connection applies in the case of a stable uniform 
liquid. 

The liquid state is described as being linearly stable if 
uj{k) < 0 for all wave numbers k and linearly unstable 
when uj{k) >0 for some wave number k. This situation 
arises for state points deep inside the parameter regime 
where the crystal is the equilibrium phase. The linear 
instability threshold is defined as the locus in the phase 
diagram where \k-k ~ ^ with uj{k = kc) = 

0, i.e., the locus where the maximum growth rate is zero. 
The location of this threshold is displayed in Fig. as 
the blue short-dashed line. 


III. STRUCTURE OF THE LIQUID 

In this section we present some typical results for the 
radial distribution function g{r), the static structure fac¬ 
tor S{k) and the dispersion relation uj{k) to illustrate the 
changes in the structure of the uniform liquid as the av¬ 
erage density po and the shoulder height parameter a are 
varied. 

In Fig. we display a series of results for fixed a = 0.8, 
as the density of the fluid is increased from zero to the 
value PqR^ = 2.7, which for a = 0.8 is the density of the 
liquid at coexistence with the crystal B phase. In the 
top panel we display the radial distribution function. In 
the limit po ^ 0 this is given by g{r) = exp[—/3U(r)], 
which exhibits a correlation hole for small r due the par¬ 
ticles seeking to avoid overlaps. However, since the re¬ 
pulsion strength for full overlap at this temperature is 
f3V{r = 0) = /3e{l + a) = 1.8, which is not that large, 
g{r « 0) is positive, reflecting the fact that there is a 
nonzero probability for particles to overlap completely, 
even at low densities. As the density po increases, the 
value of g{r ^ 0) also increases, reflecting the fact that 
particles are forced to overlap more often. Furthermore, 
oscillations develop in the tail of ^(r), at larger r. For 
particles with a hard core, we would normally ascribe this 
behavior to packing effects due to core exclusion. How¬ 
ever, in the present system this is a largely energetic 
effect: as the density is increased, the overall energy is 
lowered if some particles overlap with each other com¬ 
pletely, thereby avoiding more expensive partial overlap 
with many particles simultaneously, although the degree 
to which this occurs depends on the balance between en¬ 
ergetic and entropic effects. This behavior is also re¬ 
flected in the fact that for poR‘^ > I, g{r ^ 0) > I. This 
value of g{r ^ 0) continues to increase as the density 
is increased. For the case poR^ = 2.7 we see that g{r) 
is highly structured, with a pronounced peak at r = 0 



r/R 



kR 



FIG. 3: Correlation functions characterising the liquid phase 
for increasing density po as indicated in the key, for fixed 
a = 0.8. Top panel: the radial distribution function g{r); 
middle panel: the static structure factor S{k); bottom panel: 
the dispersion relation Cjj{k). 


indicating multiple overlaps. In fact, it is this growing 
tendency to form clusters that drives the freezing into a 
cluster-crystal when the density poR^ > 2.7. 

In the middle panel of Fig. we display the structure 
factor S{k) obtained via Eq. ([l3|) at the same density 
values. We see that as the density is increased, S{k) 
exhibits two peaks. These reflect the correlations in the 
system with two characteristic length scales, R and 
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r/R 


kR 




FIG. 4: Correlation functions characterising the liquid phase 
for a range of values of a as indicated in the key, for fixed 
poR^ = 1.2. Top panel: the radial distribution function g{r); 
middle panel: the static structure factor S{k); bottom panel: 
the dispersion relation Cjj{k). 


and so the two peaks in S{k) are (roughly) at the wave 
numbers « 27rjRs and ^ 211 jR. The fact that the peak 
at 27r/i? is larger reflects the fact that for this value 
of a the particle core repulsions dominate the repulsions 
due to the shoulder. As a result for this value of a (a = 
0.8) the system freezes to form the small lattice spacing 
crystal B phase. 

In the lower panel of Fig. |^we display the dispersion 


FIG. 5: The dispersion relation ct;(/c) for fixed p^R^ — 3.5 and 
various values of a, as indicated in the key. 



FIG. 6: The dispersion relation uj{k) at a series of points along 
a diagonal path in the phase diagram passing through the 
point at pqR^ — 2.95 and a = 1.067 at which the two modes 
kiR = 3.12 and k 2 R = 6.03 are simultaneously marginally 
unstable. Note that k 2 /ki = 1.932. 


relation uj{k) at the same series of state points. Except 
for the limiting case of low density, uj{k) also exhibits 
two peaks, reflecting the peaks in S{k). For all the re¬ 
sults displayed uj{k) < 0 for all k from which we infer 
that the liquid is in fact linearly stable at these densities. 
Indeed, at these densities the uniform liquid is the global 
minimum free energy state. However, for poR^ > 2.7, 
the global minimum corresponds to that of the hexag¬ 
onal crystal. As the density is further increased (not 
displayed), the larger k peak in uj{k) continues to grow 
in height, and when poR^ ~ 3.2 the peak growth rate 
uj{k = kc) = 0, indicating that the uniform liquid is now 
marginally unstable with respect to perturbations with 
wave number kc ^ 27r/R. The resulting linear instability 
threshold (spinodal) line is displayed as the blue short- 
dashed line in Fig. 

In Fig.we display g{r) (top), S{k) (middle) and uj{k) 
(bottom) at fixed density poR^ = 1-2, as the shoulder 

































height parameter a is varied. For small values of a, we 
see that g{r) exhibits a peak just beyond r = R, since 
this is the effective diameter of the particles. However, as 
a increases, increasing the shoulder height, this peak de¬ 
creases in height while another peak develops just beyond 
r = i?s, reflecting the growing dominance of the shoulder 
in determining the correlations in the liquid. The liquid 
with density poR^ = 1.2 and a = 1.3 is at phase coexis¬ 
tence with the crystal A phase. The behavior observed in 
g{r) is, of course, reflected in the structure factor shown 
in the middle panel of Fig. Specifically, for small a 
there is a single peak in S{k), at kR ~ 5.5. As the shoul¬ 
der height a increases, this peak moves slightly towards 
larger /c, and a second peak develops at kR ^ 3, i.e., at a 
value of k that is a little below the value 27r jRg. The lat¬ 
ter reflects the growing importance of the length scale Rs 
in the particle correlations in the liquid. As a increases 
further the peak at smaller k overtakes the larger k peak. 
The two peaks in S{k) have equal height at a = 1.067, 
irrespective of the fluid density. The locus of this point 
is displayed as the green long-dashed line in Fig. The 
lower panel of Fig. which displays the dispersion re¬ 
lation uj{k) also shows the development and growth of 
a peak at kR ^ 3. Increasing a beyond the values dis¬ 
played in this figure shows that this peak continues to 
grow in height until uj{k) > 0 for kR ~ 3, indicating 
the uniform fluid becomes linearly unstable. In Fig. we 
display the locus along which the two principal peaks in 
uj{k) are of equal height using a pink dotted line. Along 
this line the growth/decay rates for density fluctuations 
with these two wave numbers are the same. 

In Fig. we display the dispersion relation uj{k) for 
fixed poR^= 3.5 and various values of a. For the case 
a = 0.4 there is one main peak in uj{k)^ with its maximum 
close to zero, indicating that this state point is close to 
but slightly outside the linear instability threshold. As 
a increases this peak grows in height and also shifts to 
slightly larger wave numbers, as the uniform liquid be¬ 
comes linearly unstable. At the same time a second peak 
starts to develop at kR ~ 3 and becomes the dominant 
peak for a > 1.4. Since uj{k) determines the growth rate 
of density fluctuations in the unstable liquid, the figure 
reveals a transition between the fastest growing modes 
at small a to those at large a; this transition takes place 
along the pink dotted line in Fig. 

This can also be seen in Fig.[^ which displays the dis¬ 
persion relation along a diagonal path in the phase dia¬ 
gram passing through the point {pqR^^ a) = (2.95,1.067), 
corresponding to the cusp in the blue short-dashed 
marginal stability threshold line in Fig. At this point 
two modes with wavevector ratio k 2 /ki = 1.932 are 
marginally unstable. 


IV. THE CRYSTAL-LIQUID STATE 

We now turn our attention to the density profiles in the 
crystal state. As illustrated in Fig. at first sight the 



FIG. 7: (Color online) The density profile in the left panel 
of Fig. [^displayed in terms of the logarithm of the density, 
ln[p(r)^], plotted in the {x/R,y/R) plane. This representa¬ 
tion allows one to see the fine structure of the density profile 
away from the principal peaks. Note in particular the honey¬ 
comb structure surrounding each of the peaks in the crystal 
A phase on the left of the interface. 
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FIG. 8: (Color online) Top: ln[p(r)R^] in the {x/R,y/R) 
plane for a system of A = 600 particles with (a, poR^) = 
(0.8,6) confined in a square region of side L = lOR obtained 
from BD simulations (top left) and DFT (top right). The 
system forms crystal A with a density profile consisting of 
an array of peaks surrounded by a connected network within 
which the particles are free to move - this is the crystal-liquid 
state. Bottom: a snapshot from the BD simulation where each 
particle coordinate is plotted as an open circle. 


crystal A and crystal B phases appear to be standard ex¬ 
amples of hexagonally ordered cluster-crystals. However, 
closer inspection of the density profile of the crystal A 
phase reveals that this is not the case, at least at state 
points near to its coexistence with the crystal B phase 
(i.e., at smaller a values). In Fig.we display the density 
profile in the vicinity of the interface between the crys¬ 
tal A and crystal B phases shown in Fig.|^ but this time 
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FIG. 9: (Color online) (a)-(d) Plots of ln[p(r)i?^] in the crystal 
A phase in the (x/R,y/R) plane for fixed = 39 at the 
state points: (a) {a, poR‘^) = (0.75,4.1), (b) (0.9,3.8), (c) 
(1.05,3.5) and (d) (1.3, 3.1). The bottom figure shows a plot 
of the fraction of mobile particles that are in the liquid part 
of the density surrounding the density peaks. For a < 0.75 
crystal A is no longer the thermodynamic equilibrium crystal 
structure, and is replaced by crystal B. 


in terms of the logarithm of the density, ln[p(r)i7^]. This 
allows one to see the fine structure in the density profile 
in the regions of space between the main peaks of the 
hexagonal lattice. Here, we see an unbroken honeycomb¬ 
like network of density that percolates throughout the 
crystal A portion of the system, indicating that the par¬ 
ticles that contribute to this portion of the density profile 
are free to move throughout the system. To confirm the 
existence of this striking structure, we calculate, using 
both DFT and BD computer simulations, the density 
profile for a system confined within a square confining 
potential $(r) with hard walls at x = 0,1077, ^ = 0,1077 
so that ^{x,y) = 0 for (0,0) < {x,y) < (1077,1077) and 
^{x^y) = oc otherwise. The top left panel in Fig. 
shows the density profile obtained from the BD simula¬ 
tions with N = 600 particles and /7e = 1, a = 0.8, i.e., the 
average density in the box is po77^ = 6. The BD result is 
obtained simply by evolving in time the particles accord¬ 
ing to Eq. (14) and then averaging over their positions 


to calculate the density profile. The top right panel in 
Fig. shows the corresponding density profile from DFT. 
The remarkable agreement between the two confirms the 
validity of the DFT approximation for this system. The 
bottom panel in Fig. shows a snapshot showing a typ¬ 
ical configuration of the particles in the BD simulation. 
The particle positions are indicated using open circles. 
Although the majority of the particles are located on 
lattice sites, a significant minority remain mobile, with 
the particles free to move in the density lanes between 
lattice sites. The system thus consists of two dynami¬ 
cally distinct populations. This is not observed in other 
pattern-forming 2D systems, such as those in Refs. [S3- 
iD], where the dynamics of all the particles are identical. 

In Fig. bottom panel, we display the percentage of 
mobile particles in crystal A as a function of the parame¬ 
ter a, for a fixed value of the chemical potential, = 39. 
This percentage is obtained by integration over all por¬ 
tions of the density profile that are a distance 0.6577 away 
from the centre of the density peaks. Particles that con¬ 
tribute to this portion of the density are defined to be 
mobile. Figures [^a)-(d) display the logarithm of the 
density profile corresponding to the points indicated in 
the lower panel. We see that as a decreases the fraction 
of mobile particles increases from zero, reaching a value 
of over 7% at a = 0.75. We terminate the curve at this 
point because for a < 0.75 crystal A is no longer the equi¬ 
librium crystal structure. It appears that as a decreases 
below this coexistence value, the growing proportion of 
mobile particles triggers the formation of the smaller lat¬ 
tice spacing crystal B phase, whereby the mobile particles 
freeze to form the additional peaks of crystal B. 


V. THE FORMATION OF QUASICRYSTALS 

The role of resonant triads in the context of minimiz¬ 
ing a free energy with one length scale has long been 
recognized [HD ED- In particular in three dimensions 
these triads can stabilize states with icosahedral sym¬ 
metry. In two dimensions the presence of two length 
scales implies the presence of two circles of wavevectors 
in Fourier space, and resonant triads involving wavevec¬ 
tors from these two circles can also contribute to stabil¬ 
ity of quasicrystals [26], provided the interaction coef¬ 
ficients are of the correct sign. With a radius ratio of 
the two circles equal to 2cos(15°) = 1.932, equilateral 
triads, 30° triads and 150° triads involving two vectors 
from one circle and one vector from the other increase 
the number of possible triads, and so the potential con¬ 
tribution to the free energy. This configuration leads to 
dodecagonal quasicrystals; with other radius ratios, the 
situation can be yet more complicated [25|. In fact, ar¬ 
guments based on the contribution to the free energy 
from resonant triads only, important though these are, 
overlook the potential importance of higher order har¬ 
monics, whose coefficients may become arbitrarily large 
owing to the problem of small divisors that inevitably 
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appears whenever quasiperiodicity and nonlinearity oc¬ 
cur together [63]. Thus a truncation of the theory at 
cubic order, a procedure widely used in the literature, 
remains to be properly justified, although Ref. m goes 
some way towards resolving the small divisor issue. 

The mechanism identified below for stabilizing QCs in 
the present system also involves two length scales, but 
differs qualitatively from that just described (see also 
[28l [29]). In our case the system first forms the small 
length scale crystal phase. It is only when this phase 
is almost fully formed (i.e., when the dynamics is far 
into the nonlinear regime) that the longer length scale 
starts to appear, leading to the formation of the QC (see 
Figs. 10 and [IT]). Thus, what we observe is in fact a 


hitherto unseen mechanism for the formation of QCs. 


The scenario for the formation of QCs described in [28] 
requires the system (i) to be just inside the linear insta¬ 
bility line (i.e., p is restricted to a small range beyond 
Pa, the value at the linear instability line), and (ii) relies 
on the simultaneous linear growth of two distinct wave 
numbers as in [26| . In this scenario the role of the higher 
order interactions (i.e., of nonlinearity) is to stabilize the 
two length scale (QC) structures formed from the two 
linearly growing scales [28] . 


We contrast this scenario with that described here for 
a uniform liquid quenched to a region above the coexis¬ 
tence of the two crystal phases but below the pink dotted 
line in Fig. In this regime the large k peak domi¬ 
nates and small length scale density fluctuations grow 
rapidly^ig. as described by the dispersion relation 
in Fig. [T^a). In this regime the system behaves as if it 
were going to form crystal B. However, the true minimum 
of the free energy corresponds to the larger length scale 
crystal and this length scale is linearly stable (Fig.[T^a)). 
As a result, as the growing short-scale density fluctua¬ 
tions reach the nonlinear regime, the system seeks to go 
to the longer length scale structure but the smaller length 
scale imprinted from the linear growth regime leads to 
frustration. We observe this type of behavior well away 
from onset - i.e., deep inside the linear-instability thresh¬ 
old, in contrast to the scenario in [28]. In Fig. 10 we dis¬ 
play the resulting time evolution of the density profile as 
the system forms QCs. Since only one mode is unstable 
the formation of the QCs that we find can only occur via 
the nonlinear mechanism described here. 


In Fig. [^ we display DDFT results showing the for¬ 
mation of a QC structure at a = 1.067 and poR^ = 3.5 
m- The dispersion relation corresponding to this state 
point is shown in Fig. [T^b). We see that in this case the 
larger wavelength mode is no longer stable, although its 
growth rate is weak compared to that of the short wave¬ 
length mode. As a result the system first forms the pure 
small length scale crystal (see, e.g., the middle panel of 
the top row of Fig. [IT] corresponding to t* = t/rs = 2, 
where tb = /^R^/F is the Brownian timescale). How¬ 
ever, over time, starting from a grain boundary, the sys¬ 
tem evolves a QC structure much as occurs at state point 
a = 0.8, PqR^ = 3.5. In both cases this happens when 


the system is well away from the linear regime, in con¬ 
trast to the weakly nonlinear QC mechanism proposed in 
Refs. [28] [29]. Indeed, for a = 1.067 the linear instabil¬ 
ity line is at poR^ = 2.95, implying that this state point 
corresponds, like a = 0.8, poR^ = 3.5, to quite a deep 
quench. As a result, both snapshot series show that the 
linear growth regime introduces only one length scale, 
that of the small length scale crystal B phase - despite 
the presence of the weakly unstable larger length scale 
in Fig. [^ Figure also confirms that the DDFT dy¬ 
namics and the fictitious dynamics obtained from Picard 
iteration in Fig. [^ are indeed qualitatively very similar. 

As explained above, our work shows that quasicrys- 
talline structures can form even when only one of the 
two scales introduced by our choice of the potential is 
unstable; the instability forms nonlinear structures with 
this one scale only but because these do not correspond 
to the global minimum of the free energy which occurs at 
a distinct scale, the system attempts to shift the struc¬ 
ture to the thermodynamically preferred scale. This pro¬ 
cess leads to frustration that is responsible for the for¬ 
mation of the QC state. This is a qualitatively dis¬ 
tinct mechanism of QC formation from that advocated in 
Refs. [28] [29] which requires that both scales are weakly 
unstable. As a result the latter theory is only capable 
of describing QCs that have very small amplitude. In 
contrast, our quasicrystalline states are present quite far 
from the onset of instability and form from a periodic 
state via the nonlinear time-dependent process described 
above. 

The calculation of the values of a and Rs/R used above 
to home in on the parameter region where QCs might be 
observed was described in Ref. [28] as well as in earlier 
work [18]. Once the approximate parameter regime has 
been identified the details of what happens depend on 
the values of a and Rs/R. However, the QCs that we 
observe are always metastable with respect to the peri¬ 
odic crystal. By this we mean that they correspond to a 
local minimum of the free energy, but not to the global 
minimum (Fig. [Td] ). Thus the density profiles in Fig. 
are indeed possible ground states (i.e., local minima), but 
not the ground state (global minimum): the free energy 
of the state in the top panel in Fig. is slightly higher 
than that of the lower panel, but both are higher than 
that of the crystal A phase, which is the global minimum 
for this state point. As a increases beyond the range 
displayed in Fig. the QC free energy increases more 
rapidly than that of the crystal A phase - i.e., the trend 
revealed in Fig. continues and the two free energies 
do not approach one another again. In particular, the 
QC free energy is far above that of the crystal A phase 
at a = 1.067. This may also be so for the QCs obtained 
in Refs. [28] [65] [66]. In contrast, very recently [67] it 
has been shown that for the Lifshitz-Petrich free energy 
[26] . QCs are indeed the global free energy minimum for 
certain parameter values. 

In Fig. [^ we display both the QC density profiles and 
the corresponding Fourier transforms. Both exhibit 12- 
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FIG. 10: (Color online) Snapshots of ln[p(r)i?^] in the (x/R, y/R) plane obtained via Picard iteration for a = 0.8 and poR^ = 3.5, 
reve alin g the evolution towards the equilibrium state for the same state point as the results displayed in the upper panel of 
Fig. 13 The dispersion relation at this state point is displayed in Fig. 12 ^a). The panels along the top row, from left to right, 
correspond to times t = 30, 32 and 35, and along the bottom row to t = 40, 50, 200. Note that the system first forms the small 
length scale crystal (at time t ~ 30). It then tries to form the longer length scale crystal. However, due to the small length 
scale already imprinted on the system, it cannot form a perfect large length scale crystal and ends up forming a disordered 
system with domains of QC ordering. The Picard iteration used to generate these figures does not locally conserve particle 
number (although it does conserve the total density in the system - see Sec. [^, but is much faster than the full DDFT and 
gives qualitatively similar results - compare this figure with Fig. mi which is calculated with DDFT. 



FIG. 11: (Color online) Snapshots of ln[p(r)i?^] in the {x/R, v/R) plane obtained from DDFT, for a — 1.067 and p^R^ — 3.5. 
The dispersion relation at this state point is displayed in Fig. |12[ b). The panels along the top row, from left to right, correspond 
to times I/tb = t* = 1, 2 and 5, and along the bottom row to = 10, 20 and 40, where tb = fiR^/T is the Brownian timescale. 
Note that the system hrst forms a small length scale crystal (t* = 2). It then tries to form the longer length scale crystal, 
initiated from a grain boundary - see panels for t* = 5 and 10. However, because of the small length scale already imprinted 
on it, the system cannot form a perfect long length scale crystal and ends up forming a disordered system with domains of QC 
ordering - see the hnal stationary profile at =40. 
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FIG. 12: Dispersion relation at the state point a = 0.8 and 
poR^ = 3.5 (top, corresponding to Fig. [l^ and a = 1.067 
and pqR^ = 3.5 (bottom, corresponding to Fig. 11). In both 
cases, QCs form at these state points. In the upper panel 
(Fig. |10| only one mode is unstable, corresponding to the 
smaller length scale crystal B. In the lower panel (Fig. ID 
two modes are unstable, but the growth rate for the smaller 
length scale crystal B is much larger than that for the larger 
length scale crystal A. 


fold ordering. In the upper case, there is significant dis¬ 
order in the system, which is not surprising given the dy¬ 
namical mechanism we observe for QC formation. How¬ 
ever, it is possible to facilitate a more ordered final state 
by choosing (for example) a periodic domain 30 times 
larger than the shorter of the two lengthscales to allow for 
a circle of twelve vectors that are 29.98° apart and whose 
lengths differ by 0.05% [17]. Starting from an initial con¬ 
dition with these twelve modes set to a small amplitude, 
we observe that the system easily forms a ‘perfect’ exam¬ 
ple of a QC (Fig. lower panels). 

In Fig. we display density profiles obtained by tak¬ 
ing this ‘perfect’ QC and then following the solution as 
the value of Rs is changed. We find that QCs remain 
linearly stable in the Picard iteration for 1.77 < Rg/R < 
2.18. If Rs is decreased to Rg/R = 1.77, the QC solu¬ 
tion at this point becomes unstable and the Picard itera¬ 
tion leaves this solution and falls onto a crystal A profile 
(which contains defects), displayed in the left hand panel 
of Fig. If instead Rg is increased, at Rg/R = 2.03 


the Picard iteration falls off the branch of solutions corre¬ 
sponding to the ‘perfect’ QC in bottom left of Fig. [T3| onto 
a different QC branch of solutions, which is displayed in 
the middle panel of Fig. Further increasing Rg^ this 
QC then becomes linearly unstable at 77^/77 = 2.18 and 
the Picard iteration then goes to the crystal B profile 
displayed in the right hand panel of Fig. 


VI. CONCLUDING REMARKS 

In this paper we have elaborated on the results of 
Ref. [6] for a simple model soft core fluid that exhibits 
surprisingly rich phase behavior: two crystalline phases 
and a fluid phase. This stems from the fact that the pair 
potential between the particles has two length scales, 77 
and 77s, and two different energy scales, ae and (1 + a)e. 
The subtle balance of these leads to rich structuring and 
phase behavior. Pair potentials with these qualities arise 
as the effective interaction potentials between polymeric 
macromolecules. In particular, we believe that tailoring 
dendrimers with a ‘core’ plus ‘shell’ architecture should 
yield particles with effective interaction potentials akin 
to those considered here. Of course, the model system 
considered here is two-dimensional, so to observe the par¬ 
ticular behavior reported here in an experimental system, 
the particles must be confined to an interface in order to 
create an effectively two-dimensional system. The natu¬ 
ral next step to take after the work described here is to 




FIG. 13: (Golor online) Left panels: plots of ln[p(r)77^] in 
the {x/R,y/R) plane obtained from DFT for (a,po77^) = 
(0.8,3.5). Right panels: the corresponding Fourier trans¬ 
forms. The latter exhibit 12-fold symmetry, which is indica¬ 
tive of QG ordering. The top density profile is obtained from 
random initial conditions, while the lower profile was formed 
starting from an initial density profile having QG symmetry. 
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FIG. 14: (Color online) The grand potential density as a func¬ 
tion of a for fixed = 39 for the two different crystal struc¬ 
tures and also the QC solution displayed in Fig. Near 
a = 0.75 there is a point where all three have almost the 
same value of the grand potential, but the QC solution is 
never the global minimum (see inset). The crystal A phase is 
of CL type throughout the range of a shown. 

consider systems in three dimensions, where the phase 
behavior and the structures observed will be even richer. 
We are now embarking on work in this direction. 

The two most striking aspects of the present model are: 
(i) The formation of the crystal-liquid phase, having two 
dynamically distinct populations of particles, some that 
are confined to the crystal lattice sites and others that 
are mobile, residing in a honeycomb-like network around 
the main density peaks, (ii) The formation of QCs. This 
aspect is particularly interesting, because the QC struc¬ 
tures form via a mechanism that is distinct from any 
of the mechanisms that have been proposed previously. 
Namely, QC formation occurs following a deep quench 
of the uniform liquid to state points where it is unsta¬ 
ble. At these state points, the global minimum of free 
energy corresponds to the large length scale crystal A 
phase. However, in the initial linear growth regime af¬ 
ter the quench, a smaller length scale (corresponding to 
the small-length crystal B phase) grows the fastest, lead¬ 
ing to the system becoming patterned with the “wrong” 
small-wavelength density modulations. When the system 
subsequently seeks to lower its free energy and hence to 
introduce the longer length scale, it remains “stuck” with 
some ordering on the small length scale. The final equi¬ 
librium structure generally consists of a mixture of the 
two length scales and may exhibit QC ordering, i.e., the 
Fourier transform may consist of a ring of 12 peaks. The 
resulting structure is in fact a local minimum of the free 
energy, but not the global minimum. The QCs formed 
via such a mechanism are, unsurprisingly, generally dis¬ 
ordered, containing a mixture of domains with 12-fold or¬ 
dering and domains of hexagonal ordering, corresponding 
to one or other of the two hexagonal crystal structures. 

The results presented here are for just one tempera¬ 
ture. However, the important quantities for determin¬ 


ing the phase behavior of the model are the dimension¬ 
less quantities ksTje^ a and Rs/R. Varying these de¬ 
termines the location in the phase diagram of the lin¬ 
ear instability threshold, the point where both length 
scales are marginally unstable and the ratio k 2 /k\. In 
the limit a = 0, increasing ksT/e shifts the linear in¬ 
stability threshold to higher density p [55]. Varying the 
temperature by a modest amount should leave the behav¬ 
ior of the present system qualitatively unchanged, merely 
shifting the regions where the crystalline phases occur to 
higher densities. 

It is worth connecting the present work with related 
work [55l [68] [69] on the freezing of binary mixtures of 
particles. The mixtures considered in Refs. [55] [68] also 
possess two length scales, owing to the fact that they 
are a binary mixture of soft particles of different sizes, 
and form multiple structures when a solidification front 
advances into an unstable uniform liquid. For a deep 
enough quench, such a front deposits behind it density 
modulations that are also of the “wrong” wavelength, 
thereby frustrating the formation of a well-ordered “cor¬ 
rect” wavelength equilibrium crystal. In particular, the 
final equilibrium structures also contain a high degree 
of disorder. In this case, the selection of the “wrong” 
wavelength is due to the dynamical nature of the length 
scale selection problem via an advancing front: the se¬ 
lected wavelength depends only on the linearization 0, 
whereas the global minimum free energy crystal structure 
is determined by the full DFT, which is highly nonlinear. 
This situation differs from the QC formation observed 
in the present work, yet there are similarities: both sys¬ 
tems undergo a linear process that generates modula¬ 
tions with a length scale that does not correspond to the 
length scale of the equilibrium structure, which is deter¬ 
mined by nonlinear processes. This naturally leads to 
an unanswered question: what happens when a solidifi¬ 
cation front advances in the present system? The front 
motion will generate a particular length scale, the linear 
growth of any local density modulations will produce a 
slightly different length scale while nonlinear interactions 
will seek to generate a third length scale. We anticipate 
that the interplay of such processes will inevitably lead 
to disordered structures. 
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